#include "similar.h"

int32_t
run_lcs_similar(char* s1, uint32_t sz1, char* s2, uint32_t sz2, char **lcs)
{
  int32_t **M;
  int32_t *lcs_a;
  int32_t diag_pos, len_max, len_lcs, pos_1, pos_2, 
    lower_2, upper_2, pos_lcs, i;

  if (sz1 > sz2) {
    uint32_t sztemp = sz2;
    char *stemp = s2;
    s2 = s1;
    sz2 = sz1;
    s1 = stemp;
    sz1 = sztemp;
  }

  M = (int32_t **) malloc(sizeof(int32_t *) * (sz1 + 1));

  for (int32_t j = 0; j < (int32_t)(sz1 + 1); j++) {
      M[j] = (int32_t *) malloc(sizeof(int32_t) * (sz1 + 1 - j));
      // NOTE: Think we can get away with this?
      //      memset(M[j], -1, sizeof(int32_t) * (sz1 + 1 - j));
  }

  diag_pos = sz1 - 1;
  len_max = 0;
  pos_1 = upper_2 = -1;

  // Calculate M
  while (len_max < diag_pos) {
    pos_1 = diag_pos;
    len_lcs = 0;
    upper_2 = sz2;
    
    // Calculate M[i][j] on one diagonal
    while (pos_1 != -1 && upper_2 != -1) {
      if (diag_pos == (int32_t)(sz1 - 1) || len_lcs > len_max) {
	M[len_lcs][pos_1 + 1] = -1;
      }

      lower_2 = MAX(0, M[len_lcs][pos_1 + 1]);
      pos_2 = upper_2 - 1;

      // Calculate one M[i][j] according to Lemma 3
      while (pos_2 >= lower_2 && s1[pos_1] != s2[pos_2]) {
	pos_2--;
      }

      if (pos_2 >= lower_2) {
	upper_2 = pos_2;
      }
      else {
	upper_2 = M[len_lcs][pos_1 + 1];
      }

      M[len_lcs][pos_1] = upper_2;

      if (upper_2 == -1) {
	len_lcs--;
      }

      len_max = MAX(len_max, len_lcs);
      len_lcs++;
      pos_1--;
    }

    diag_pos--;
  }

  lcs_a = (int32_t *) malloc(sizeof(int32_t) * (sz1 + 1));
  memset(lcs_a, -1, sizeof(int32_t) * (sz1 + 1));

  // Recover LCS
  if (upper_2 == -1) {
    pos_1 += 2;
  }
  else {
    pos_1++;
  }

  pos_lcs = len_max;

  while (pos_lcs > -1) {
    while (M[pos_lcs][pos_1] == M[pos_lcs][pos_1 + 1]) {
      pos_1++;
    }

    lcs_a[pos_lcs] = pos_1;
    pos_1++;
    pos_lcs--;
  }

  for (i = 0; i < (int32_t)(sz1 + 1); i++) {
    free(M[i]);
  }

  free(M);

  len_lcs = 0;
  for (i = 0; i < (int32_t)(sz1 + 1); i++) {
    if (lcs_a[i] == -1) break;
    len_lcs++;
  }

  *lcs = (char *) malloc(sizeof(char) * len_lcs);
  
  for (i = 0; i < len_lcs; i++) {
    (*lcs)[len_lcs - 1 - i] = s1[lcs_a[i]];
  }

  free(lcs_a);

  return len_lcs;
}
